/*
This Stata do file replicates the main results from Kinne, Brandon J. 2012. "Multilateral
Trade and Militarized Conflict," in Journal of Politics. The estimates used in Figure 3
of the paper are generated by this do file but are not plotted within Stata. Email the
author (brandon.kinne@utdallas.edu) for an R script to plot these estimates.

All results in the paper were generated using this script in Stata/MP 11.2 on a 64-bit
Intel Mac running OS X 10.6. The script was also checked on 64-bit Ubuntu.
*/

/* Open the do file from the directory in which Kinne_JOP.dta is located */

use Kinne_JOP, clear

tsset ccode year, yearly

summ

***********
* TABLE 2 *
***********

/* Set up the estimation */
local base = "cent50 lnopen polity lngdppc lncinc contigs concap" /* List of baseline controls, plus centrality */
capture drop sample
quietly xtgee ftmids `base', family(nbinomial) link(nbinomial) corr(ar 1) vce(robust)
quietly gen sample = e(sample) /* Record the GEE sample so same N can be used for logit */

/* Baseline Model */
quietly nbreg ftmids `base' if sample==1
local a = e(alpha) /* Capture alpha */
xtgee ftmids `base', family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

/* With lag of DV */
quietly nbreg ftmids `base' totmids if sample==1
local a = e(alpha)
xtgee ftmids `base' totmids, family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

/* Fixed effects negative binomial */
xtnbreg ftmids `base', fe

/* ART model */
quietly nbreg ftmids cent50 lnopen lncinc if sample==1	/* Would probably be better here to set alpha using larger sample... */
local a = e(alpha)
xtgee ftmids cent50 lnopen lncinc, family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

/* Logit models with binary DV */
logit mid `base' pyrs _spline1 _spline2 _spline3 if sample==1, vce(cluster ccode)

/* Fixed effect logit with binary DV */
xtlogit mid `base' pyrs _spline1 _spline2 _spline3, fe noskip

***********
* TABLE 3 *
***********

/* Trade power in place of centrality */
quietly nbreg ftmids lntpower lnopen polity lngdppc lncinc contigs concap if sample==1
local a = e(alpha)
xtgee ftmids lntpower lnopen polity lngdppc lncinc contigs concap, family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

/* Then both together */
quietly nbreg ftmids cent50 lntpower lnopen polity lngdppc lncinc contigs concap if sample==1
local a = e(alpha)
xtgee ftmids cent50 lntpower lnopen polity lngdppc lncinc contigs concap, family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

/* Now residualization */
reg cent50 lntpower
predict r_cent, r
label var r_cent "Centrality residual"

quietly nbreg ftmids r_cent lntpower lnopen polity lngdppc lncinc contigs concap if sample==1
local a = e(alpha)
xtgee ftmids r_cent lntpower lnopen polity lngdppc lncinc contigs concap, family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

reg lntpower cent50
predict r_tpow, r
label var r_tpow "TradePower residual"

quietly nbreg ftmids cent50 r_tpow lnopen polity lngdppc lncinc contigs concap if sample==1
local a = e(alpha)
xtgee ftmids cent50 r_tpow lnopen polity lngdppc lncinc contigs concap, family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)

*****************
* FOR R FIGURES *
*****************

/* First scale cent50 and lnopen between zero and one, to make graphing easier */
quietly summ cent50
replace cent50 = cent50 / r(max)

quietly summ lnopen
replace lnopen = lnopen / r(max)


gen yrgrp5 = .	/* Label groups for regressions on five-year subsamples */
forvalues i = 1954(5)1989 {	/* GEE model doesn't converge for 1950-1954 period, so drop that one */
	local j = `i' + 6
	replace yrgrp5 = `i'+1 if year>`i' & year<`j'
}
replace yrgrp5 = 1995 if year>1994

gen indx = .
egen drops = fill(1955 1960 1965 1970 1975 1980 1985 1990 1995)
forvalues i = 1(1)9 {
	replace indx = drops[`i'] in `i'
}
drop drops

local base = "cent50 lnopen lncinc" /* Use ART model for five-year subsamples */

/* GEE estimates */
gen mngee = . /* Estimate for centrality */
gen segee = . /* SE for centrality */
gen omngee = . /* Estimate for openness */
gen osegee = . /* SE for openness */
forvalues i = 1955(5)1995 {
	quietly nbreg ftmids `base' if yrgrp5==`i'
	local a = e(alpha)
	capture xtgee ftmids `base' if yrgrp5==`i', family(nbinomial `a') link(nbinomial) corr(ar 1) vce(robust)
	replace mngee = _b[cent50] if indx==`i'
	replace segee = _se[cent50] if indx==`i'
	replace omngee = _b[lnopen] if indx==`i'
	replace osegee = _se[lnopen] if indx==`i'
}

/* Logit estimates */
gen mnlog = .
gen selog = .
gen omnlog = .
gen oselog = .
forvalues i = 1955(5)1995 {
	quietly logit mid `base' pyrs if yrgrp5==`i', vce(cluster ccode)
	replace mnlog = _b[cent50] if indx==`i'
	replace selog = _se[cent50] if indx==`i'
	replace omnlog = _b[lnopen] if indx==`i'
	replace oselog = _se[lnopen] if indx==`i'
}

/* Predicted probabilities from the Logit model */
gen prmid = .
gen prsd = .
gen oprmid = .
gen oprsd = .
/* First do centrality */
forvalues i = 1955(5)1995 {
	estsimp logit mid `base' pyrs if yrgrp5==`i', vce(cluster ccode)
	setx mean
	quietly summ cent50 if yrgrp5==`i'
	setx cent50 r(mean)-r(sd)
	quietly simqi, prval(1) genpr(pr1)
	quietly summ cent50 if yrgrp5==`i'
	setx cent50 r(mean)+r(sd)
	quietly simqi, prval(1) genpr(pr2)
	gen delta = ((pr2-pr1)/pr1)*100
	summ delta
	replace prmid = r(mean) if indx==`i'
	replace prsd = r(sd) if indx==`i'
	drop pr1 pr2 delta
	drop b*
}
/* Then do openness */
forvalues i = 1955(5)1995 {
	estsimp logit mid `base' pyrs if yrgrp5==`i', vce(cluster ccode)
	setx mean
	quietly summ lnopen if yrgrp5==`i'
	setx lnopen r(mean)-r(sd)
	quietly simqi, prval(1) genpr(pr1)
	quietly summ lnopen if yrgrp5==`i'
	setx lnopen r(mean)+r(sd)
	quietly simqi, prval(1) genpr(pr2)
	gen delta = ((pr2-pr1)/pr1)*100
	summ delta
	replace oprmid = r(mean) if indx==`i'
	replace oprsd = r(sd) if indx==`i'
	drop pr1 pr2 delta
	drop b*
}

/* Send these figures to an R script or wherever you like to do graphics */
list mngee segee omngee osegee mnlog selog omnlog oselog prmid prsd oprmid oprsd in 1/9
